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ABSTRACT 

We describe simple useful toy models for key processes of galaxy formation in its 
most active phase, at z > 1, and test the approximate expressions against the typi- 
cal behaviour in a suite of high-resolution hydro-cosmological simulations of massive 
galaxies at z — 4 — 1. Wc address in particular the evolution of (a) the total mass 
inflow rate from the cosmic web into galactic haloes based on the EPS approximation, 
(b) the penetration of baryonic streams into the inner galaxy, (c) the disc size, (d) the 
implied steady-state gas content and star-formation rate (SFR) in the galaxy subject 
to mass conservation and a universal star-formation law, (e) the inflow rate within 
the disc to a central bulge and black hole as derived using energy conservation and 
self-regulated Q ^ I violent disc instability (VDI), and (f) the implied steady state in 
the disc and bulge. The toy models provide useful approximations for the behaviour 
of the simulated galaxies. We find that (a) the inflow rate is proportional to mass and 
to (1 -I- z)^/^, (b) the penetration to the inner halo is ^ 50% at z = 4 — 2, (c) the disc 
radius is ^ 5% of the virial radius, (d) the galaxies reach a steady state with the SFR 
following the accretion rate into the galaxy, (e) there is an intense gas inflow through 
the disc, comparable to the SFR, following the predictions of VDI, and (f) the galaxies 
approach a steady state with the bulge mass comparable to the disc mass, where the 
draining of gas by SFR, outflows and disc inflows is replenished by fresh accretion. 
Given the agreement with simulations, these toy models are useful for understanding 
the complex phenomena in simple terms and for back-of-the-envelope predictions. 

Key words: cosmology — dark matter — galaxies: evolution — galaxies: formation 
— galaxies: haloes 



1 INTRODUCTION 

A scientiflc understanding of a physical phenomenon is 
often materialized through a simple, idealized, analytic 
model, which we call a toy model. A toy model is de- 
rived from first principles in an attempt to capture the 
key physical elements of the process and to provide a 
useful approximation to the complex behaviour seen in 
numerical simulations or in observations. The simple 
scaling relations from such a toy model allow back-of- 
the-envelope calculations that involve several elements 
or a sequence of events, which may help to develop an 
understanding for a whole scenario, and hopefully lead 
to explanations of observed phenomena or to new the- 
oretical predictions. 

A different tool in the study of galaxy formation 
is semi-analytic modeling (SAM , e.g.,lSomerville et al.l 
[2008t iGuo et al.l [20ll iBensonI 120]^ . In SAMs, the 



many complex physical processes associated with gas, 
stars, black holes and radiation are modeled in par- 
allel by recipes with free parameters, incorporated in 
retrospect in dark-matter halo merger trees that were 
obtained from dissipationless cosmological N-body sim- 
ulations o r the Extended P ress-Schechter approxima- 
tion (EPS. lBond et al.lll99lT ). The toy models addressed 
here are much simpler than the sophisticated SAMs. 
The toy models are typically based on simple analytic 
arguments using crude approximations, and they tend 
to focus on one process at a time. These toy models 
could serve as a basis for physical recipes to be incor- 
porated in SAMs. SAM recipes have been confronted 
with the more elaborate hydro-cosmological simula- 
tions, where gas processes are explicitly simulated (e.g. 
Benson et al.ll200ll:ICattaneo et~aLll2007l: fNeistein etall 



20121 : iHirschmann et al.ll2012t ). but this is commonly 
done in the context of a full SAM. Here we compare 
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each toy model on its own to the corresponding feature 
of the simulation results. 

We address toy models for the formation and evolu- 
tion of galaxies in their most active phase, at redshifts 
z > 1, where the Einstein-deSitter (EdS) cosmological 
model serves as a useful approximation. We focus on 
central galaxies in their haloes, the main progenitors of 
central galaxies at z 1. In particular, we address sev- 
eral key processes as follows, (a) We verify the inflow 
rate of dark matter and baryons from the cosmic web 
into distinct haloes based on the EPS approximation, 
(b) We evaluate the penetration of baryonic streams 
through the hot gas in the haloes into the inner disc 
galaxies, (c) We crudely estimate the disc radius as a 
function of the virial radius compared to the standard 
model where angular momentum is conserved, (d) We 
study the implied gas content and star-formation rate 
(SFR) in the galaxy subject to mass conservation (the 
"bathtub" model) and a universal SFR law. (e) We ad- 
dress the inflow rate within the disc to a central bulge 
as derived using energy conservation and Toomre Vi- 
olent Disc Instability (VDI). (f) We study the implied 
steady state in the disc as it is fed by cold streams and 
is drained by star formation, outflows, and inflow to the 
bulge. 

To test the validity of the toy models, we com- 
pare their predictions with the results obtained from 
a suite of 27 Adaptive Mesh Refinement (AMR) hydro- 
cosmological simulations that zoom in on individual 
galaxies with a resolution of ^ 50 pc. These are state- 
of-the-art simulations that resolve the key processes ad- 
dressed. The feedback in these simulations is limited to 
the thermal effects of stellar winds and supernovae, not 
addressing radiative feedback and AGN feedback, thus 
not allowing an immediate study of the effects of very 
strong outflows and the associated severe suppression 
of SFR. 

The outline of this paper is as follows. In ^21 we de- 
scribe the simulations. In S|3]wc address the mass inflow 
rate into haloes. In 21 we study the penetration of cold 
streams into the central galaxy at the inner halo. In fj5] 
we refer to the disc size as derived from conservation 
of specific angular momentum. In fJHlwc investigate the 
steady-state of the gas and stellar content in which the 
SFR follows the accretion rate. In 33 we address the 
VDI-driven inflow within the disc into the bulge. In ^J5] 
we refer to the resultant steady state in the disc. In 
we summarize our conclusions. 



2 SIMULATIONS 

We use zoom-in hydro cosmological simulations 
of 27 moderately massive galaxies with an AMR 
maximum resolution 35 — 70 pc, all evolved from 
high redshift to 2 = 2, many to 2; 1.5 and 
some reaching z — 1, as marked by agn in Ta- 
ble [TJ They ut ilize the Adaptive Refinemen t Tree 
(ART) code (iKravtsoy. Klvpin fc Khokhlovl 119971: 
ICeverino fc KlvpinI 1200^ " which accurately follows 
the evolution of a gravitating N-body system and 



the Eulerian gas dynamics using an adaptive mesh. 
Beyond gravity and hydrodynamics, the code incor- 
porates at the subgrid level many of the physical 
processes relevant for galaxy formation. They include 
gas cooling by atomic hydrogen and helium as well 
as by metals and molecular hydrogen, photoioniza- 
tion heating by the UV background with partial 
self-shielding, star formation, stellar mass loss, metal 
enrichment of the ISM, and feedback from stellar 
winds and supernovae, im plemented as local injec 
tion of thermal energy JCeverino fc Klvpin 200 
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2.1 Subgrid Physics 

A few relevant details concerning the subgrid physics 
are as follows. Cooling and heating rates are tabu- 
lated for a given gas density, temperature, metallic- 
ity and UV backgro und based on the CLOUDY code 
(|Ferland et al.lll998D . assuming a slab of thickness 1 
kpc. A uni form UV background ba sed on the redshift- 
dependent iHaardt fc Madaul (|1996[ ) model is assumed, 
except at gas densities higher than 0.1 cm~'^, where a 
substantially suppressed UV background is used (5.9 x 
10^^ergs~^cm~^Hz~^) in order to mimic the partial 
self-shielding of dense gas. This allows the dense gas 
to cool down to temperatures of ~ 300K. The as- 
sumed equation of state is that of an ideal mono-atomic 
gas. Artificial fragmentation on the cell size is pre- 
vented by introducing a pressure floor, which ensures 
that the Jeans scale is resolved by at least 7 cells (see 
ICeverino. Dekel fc Bournaudll2010[ ). 

Star formation is assumed to occur at cell densi- 
ties above a threshold of 1 cm^'^ and at temperatures 
below lO^K. More than 90% of the stars form at tem- 
peratures well below lO^K, and more than half the stars 
form at 300 K in cells where the gas density is higher 
than 10 cm"'^. The code implements a stochastic star- 
formation model that yields a star-formation efficiency 
per free-fall time of 5%. At the given resolution, this 
efficiency roughly mimics b y construction the empir- 
ical K ennicutt- Schmidt law ()Schmidtlll959t iKennicuttj 
Il998l) . When extracting the resultant approximate SFR, 
we consider the mass in stars born in a given timestep 
divided by the elapsed time. The code incorporates a 
thermal stellar feedback model, in which the combined 
energy from stellar winds and supernova explosions is 
released as a constant heating rate over 40 Myr follow- 
ing star formation, the typical age of the lightest star 
that explodes as a type-II supernova. The heating rate 
due to feedback may or may not overcome the cool- 
ing rate, depend ing on the gas conditi ons in the star- 
formi ng regions (jDekel fc Si lk 1986; Ce verino fc KlvpinI 
I2OO9I) . We note that no shutdown of cooling or any other 
artificial mechanism for boosting the feedback effects is 
implemented in these simulations. We also include the 
effect of runaway stars by assigning a velocity kick of 
^ lOkms"^ to 30% of the newly formed stellar parti- 
cles. The code also includes the later effects of type-la 
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Table 1. The suite of simulated galaxies. Quoted at a = 0.33 (z = 2) are the virial radius Rv and within it the total mass Mv, stellar 
mass J\/star and gas mass Mg. The snapshots used are from a = 0.2 (z = 4) to o = afin (^Sn = a,'^^ — 1 = 2 — 1), and the virial mass at 
afi„ is quoted. 




Figure 1. Growth of virial mass Mv as a function of expansion factor a = (1 + z) ^ for the individual simulated galaxies. Each galaxy 
has been evolved to a final redshift as specified in Table [T] The simulations span about a decade in mass, from Mv = 0.16 X 10^2 to 
1.9 X IO^^Mq at z = 2. The solid black curves refer to the toy-model prediction, eq. l(9]l, normalized to Mv = W^^Mq at z = 2. The 
overall growth pattern is reproduced by the toy model quite well, though the individual galaxies can grow in different rates at different 
times in their histories. 



supernova and stellar mass loss, and it follows the metal 
enrichment of the ISM. 

2.2 Selected Haloes 

The initial conditions for the high-resolution, zoom-in, 
hydrodynamical simulations that are used in this pa- 



per were based on dark-matter haloes that were drawn 
from dissipationless N-body simulations at lower reso- 
lution in three large comoving cosmological boxes. The 
assumed cosmology is the standard ACDM model with 
the WMAP5 values of the cosmological parameters, 
namely ^ = 0-27, 17a = 0.73, = 0.045, h = 0.7 and 
(78 = 0.82 (jKomatsu et al.l [20091 ) . Distinct haloes were 
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identified about density peaks, and their virial radius 
and mass, i?v and Mv, were measured from tlic radial 
mass profile such that the mean overdensity within the 
sphere of radius i?v equals A(a) ( ^A21 eq. IA8|) . We do 
not explicitly address the evolution of subhaloes in this 
work, but rather consider them as part of the accre- 
tion onto the the central galaxy of their host distinct 
halo. The haloes for re-simulation were selected to have 
a given virial mass at z = 1 (or z = in a few cases) , as 
specified in Table [T] The only other selection criterion 
was that the haloes show no ongoing major merger at 
z = 1. This eliminates from the sample less than 10% 
of the haloes which tend to be in a dense environment 
at 2 '--^ 1, and it induces only a minor selection effect at 
higher redshifts. The target virial masses at z — I were 
typically selected to be Afv (1-5 - 3) x IO^^Mq, with 
most of them intended to end up as (3 — 6) x IO^^Mq to- 
day if left in isolation, namely somewhat more massive 
than the Milky Way. However, the actual mass range 
of these haloes at z = is broader, with some of the 
haloes destined to merge into more massive haloes host- 
ing groups or clusters. 

2.3 Zoom in 

The initial conditions corresponding to each of the se- 
lected haloes were filled with gas and refined to a much 
higher resolution on an adaptive mesh within a zoom-in 
Lagrangian volume that encompasses the mass within 
twice the virial radius at z = 1, which is roughly a 
sphere of comoving radius 1 Mpc. This was embed- 
ded in a comoving cosmological box of side 20, 40 or 
80 h~^Mpc. Each galaxy has been evolved with the full 
hydro ART and subgrid physics on an adaptive comov- 
ing mesh refined in the dense regions to cells of mini- 
mum size between 35-70 pc in physical units at all times. 
This maximum resolution is valid in particular through- 
out the cold discs and dense clumps, allowing cooling 
to ^ 300K and gas densities of ~ lO'^cm"'^. The dark- 
matter particle mass is 6.7 x IO^Mq, and the particles 
representing groups of stars have a minimum mass of 
IO^Mq. 

The virial properties of all 27 galaxies in our sam- 
ple are listed in Table [TJ This includes the virial ra- 
dius and total virial mass at z = 2, the stellar mass 
and gas mass within the virial radius at that time, the 
latest time of analysis for each galaxy in terms of the 
expansion factor, afin, and the virial mass at that time. 
The g alaxies MWOl-03 were the bas i s for the study of 
VDI in lCeverino. Dekel fc BournaudI (|2010D . The galax- 
ies MW04 and SFGl w ere added to th e stud y of giant 
clumps in VDI discs bv lCeverino et al.l (|2012[ ). 

2.4 Analysis 

We start the analysis presented below at the cosmo- 
logical time corresponding to expansion factor a = 0.2 
(redshift z = 4). At earlier times, the fixed resolution 
scale typically corresponds to a larger fraction of the 
galaxy size, which may bias some of the quantities that 
we wish to study here. All 27 galaxies reach a = 0.33 



(z = 2), 17 galaxies reach a = 0.4 (z = 1.5), and only 7 
galaxies have been run all the way to a — 0.5 (z = 1). 
This gradual degradation of the sample after z = 2 has 
been taken into account in our analysis. 

The output of each simulation is provided at out- 
put times separated by a constant interval in a, Aa, 
which for most galaxies is Aa = 0.01. For two galaxies 
(SFG8-9) the timestep is twice as small, Aa = 0.005, 
and for four galaxies (MWOl-4) the timestep is twice as 
large, Aa = 0.02. For every galaxy we analyze the data 
averaged over timesteps of Aa = 0.02. 

The quantities addressed, for example, are mass M, 
inflow rate M, and specific inflow rate M/M. For com- 
parison with a toy model at a given time step, we stack 
the data from all the available galaxies, sometimes af- 
ter proper scaling motivated by the model in order to 
minimize possible mass dependence. For each quantity, 
wc compute and show the median and the linear aver- 
age. Because of the lognormal nature of the distribution 
(see below, tj3.3.21 Fig. [S]), we also compute the aver- 
age of the log within the 90% percentiles. The tails are 
eliminated here because M could be negative in a few 
rare cases due, for example, to a fly-by satellite caught 
while crossing the virial radius outwards. The exclusion 
of the outliers also helps reducing fluctuations between 
time steps. The scatter is typically marked in our figures 
by shaded areas between the upper and lower 68% per- 
centiles and 90% percentiles. We also show as error-bars 
an estimate for the error of the mean, a/y/N, where a 
is the standard deviation over the N galaxies averaged 
over at the given time step . 

The statistical analysis considers all the galaxies of 
the sample in the time range a = 0.2 — 0.33 (z = 4 — 2), 
but only the gradually degrading subsample between 
a = 0.33 (z = 2) and a = 0.5 (z = 1) according to 
the value of afin in Table [TJ This translates to larger 
uncertainties in the results after z = 2. We therefore 
focus most of our attention on the simulations at z ^ 2. 
When addressing detailed disc properties (e.g., in m.2\ 
and iJS]), we limit the main analysis to the range z = 
3-2. 

2.5 The Sample of Galaxies 

Figure [T] shows the curves of virial mass growth in time 
for the different galaxies listed in Table [1] At z = 2, 
the sample spans roughly an order of magnitude in halo 
mass, from Mv = 0.16 x 10" to 1.9x IO^^Mq, with virial 
radii from 70 to 135 kpc. The simulated haloes typi- 
cally contain within the virial radius a baryon fraction 
of 0.14, slightly smaller than the cosmic value, of which 
the stellar fraction is typically 0.10. This suite of galax- 
ies spans the mass range of t ypical observed massive 
star-f orming galaxies at z ^ 2 (jForster Schreiber et al.l 
l2009j ). and otherwise no major selection criteria was 
i mposed on their properties at z ^ 2 — 4. As shown 
m ICevermo. Dekel fc BournaudI (|2010D . our first simu- 
lated galaxies, MWOl-03, are consistent with the ob- 
served scattered scaling relations of z '--^ 2 galaxies, in- 
cluding the relation betw een SFR and stellar mass and 
the TuUy-Fisher relation (jForster Schreiber et aTll2009l : 
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ICresci et al.l l2009l but see a discussion of deviations 
below). We can therefore assume that this is approx- 
imately a fair sample of galaxies in the relevant mass 
and redshift range, excluding galaxies that will be in 
rather dense environments at z ~ 1. 

The mass growth curves shown in Fig. [T] are to be 
visually compared to the toy model prediction discussed 
later in eq. ©, which is arbitrarily normalized in this 
figure to = lO^^Mg at z = 2. The overall growth 
pattern is similar to the toy model prediction, though 
the individual galaxies can grow in different rates at 
different times in their histories. We note that in some 
cases galaxies tend to keep their rank order with respect 
to the other galaxies in terms of virial mass, but in other 
cases galaxies change their rank order drastically during 
the evolution (e.g., MW03, left panel, sohd green line), 
commonly following major mergers in their histories. 
The simulated growth rates are compared to the toy 
model in more detail in Figs. [2] to H 



2.6 Limitations of the Current Sample 

These simulations are state-of-the-art in terms of the 
high-resolution AMR hydrodynamics and the treatment 
of key physical processes at the subgrid level. In par- 
ticular, they properly trace the cosmological streams 
that feed galaxies at high redshift, including merg- 
ers and smooth flows, and they resolve the violent 
disc instability that gove rns the high-z disc evolution 
and t he bulge format i on (ICeverino. Dekel fc BournaudI 
l2010t ICeverino etall |2012[ ). When tracing the im- 
portant small-scale processes involved in galaxy evolu- 
tion, AMR codes offer a distinct advantage over SPH 
codes, e.g., because SPH does not accurately model 
processes such as shocks and hydrodynamical instabil- 
ities (e.g. lAgertz et al.l 120071 : IScannapieco et all 120121 : 
iBauer fc Springell |2012[ ). The AMR codes seem to be 
comparable in their capa bilities to new codes us ing a 
moving unstructured grid ([Bauer fc Springeill2012[ ) , but 
the latter are currently limited to signif icantly lower 
resolu tion, on the order of ~ Ikpc (e.g., iNelson et all 
l2013f ). which is not sufficient for tracing the key pro- 
cesses of disc instabilities or the evolution of streams. 

However, like other simulations, the current 
simulations do not yet treat the star formation 
and feedback processes with sufficient accuracy. 
For example, the code assumes a somewhat high 
SFR efficiency per free-fall time, it does not fol- 
low in detail the formation of molecules and the 
effect of metallicity on SFR (jKrumholz fc Dekell 
l2012f ). and it does not explicitly include radiative 
stellar feedback (iMurrav. Quataert fc Thompson 



201C; 



20121: 



IKrumholz fc Dekell 
IDekel fc K rumholj mi^ 
(ISilk fc Re es 1991 



20101 iHopkins et al 



or 



_ AGN feed- 

back (ISilk fc Re es 19981: IHopkins et all 120061: 
iBooth fc Schavell2009l : ICattaneo et al.ll2009D . Therefore, 
the early SFR is overestimated, while the suppression 
of SFR in small galaxies is underestimated, resulting 
in excessive early star formation prior to z ~ 3, by a 
factor of order 2. As a result, the typical gas fraction 
and SFR at z ^ 2 arc lower by a factor of ^ 2 than 



the average observed values in star-forming galaxies 
(Ceverino, Dekel fc BournaudI l2010l : iDaddi et al.ll2010l: 
pTacconi et allfeOlOi ). 

Furthermore, the simulated galactic mass outflow 
rate is only a fraction of the SFR, where the mass 
loading factor ranges from zero to unity with an av- 
erage 7] ~ 0.3 at 0.5i?v, not reproducing some of 
the observed strong outflows with mass l oading fac- 
tors of order unity and above jS teidel et al.l 120101 : 
iGenzel et aLll201ll : lDekel fc Krumholzii2013l) . This leads 
to a stellar fraction of ~ 0.1 within the virial ra- 
dius, a factor of ^ 2 — 3 higher than the observation- 
ally indicated value (e.g. Perez-Gonzalez et al.l 120081 : 



iBehroozi. Wechsler fc Conrovl 12013^ These inaccura- 
cies in the SFR, feedback and outflows introduce a lim- 
itation on the generality of our testing. Nevertheless, 
this situation has the advantage that one can test the 
toy models for cosmological accretion and the response 
of the galaxy to it without the extra complication of 
very strong feedback effects. More accurate recipes for 
star formation and feedback are being incorporated into 
simulations that we and others are now running, and 
they will be improved further in the future. This will 
enable the next generation comparison to toy models. 

One should note that the sample of simulated galax- 
ies is not a fair sample of halo or galaxy mass. The halo 
masses arc limited by selection to a rather narrow mass 
range at the target redshift z ~ 1, and therefore to lim- 
ited mass ranges at earlier redshifts, and the mass dis- 
tribution at any redshift does not follow the ACDM halo 
mass function. We do not attempt to follow the proper- 
ties of a population of galaxies as it evolves, partly be- 
cause of the mass variation at z = 1, and partly because 
of the degrading of our sample at late times. Instead, 
we treat each outputcd snapshot independently of its 
history, assuming that it represents an (almost) arbi- 
trary galaxy of an instantaneous mass M at z. Since the 
galaxies are of different masses at any given redshift, we 
evaluate "average" properties for the simulated galaxies 
focusing on quantities that arc only weakly dependent 
on mass within the spanned mass range, such as the 
specific accretion rate. For quantities that do have a 
significant mass dependence, such as the accretion rate 
and the cumulative mass growth, we attempt to scale 
out the mass dependence using the toy model itself prior 
to stacking the galaxies of different masses. 

It should be emphasized that we only follow the 
growth of the central galaxies in the main-progenitor 
haloes of the final haloes selected at z = 1. The other 
progenitors, satellite and merging haloes (or galaxies), 
arc all considered here as part of the mass inflow onto 
the main progenitor. The growth of the non-main- 
progenitor galaxies in the zoom-in simulations is not 
studied here, partly because they represent only 10-20% 
of the population of galaxies of the same mass, which 
mostly consists of central galaxies, not satellites of more 
massive galaxies. We do not attempt here an extension 
of the zoom-in simulations to smaller central galaxies 
because the accuracy at the given resolution becomes 
limited, and because the proper treatment of feedback 
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effects in small galaxies, where they dominate the evo- 
lution, is more demanding. 



3 COSMOLOGICAL ACCRETION RATE 

3.1 Halos in the Einstein-deSitter Regime 

We consider the ACDM cosmology in the Einstein- 
deSitter (EdS) regime. This is a useful approximation at 
z > 1, and it becomes more and more accurate at higher 
redshifts. In our toy modeling, in the EdS regime, we 
adopt the following simple relation between the expan- 
sion factor a and the Hubble time t. 



a = (1 + zY 



2/3 



(1) 



ti = -r!,-i/2i/-i^l7.5Gyr. (2) 

A comparison to the accurate expressions in the ap- 
pendix, reveals that this approximation overesti- 
mates t a.t z ^ 1 hy only 5%, and by 1.6% at z = 2 
(but by 28% at z = 0). The mean mass density in the 
Universe is 



Pu = Poa 



2.5 X 10"^°gcm-3 . 



(3) 



Inspired by the spherical-collapse model and the 
virial theorem, a halo is defined as a sphere about a den- 
sity peak that in the EdS regime encompasses a mean 
overdensity of A ~ 200 above the cosmological back- 
ground (see eq. (|A8p for a more accurate expression, 
where A ~ 207, 187 and 178 at z = 1,2 and z > 1). 
The relations between virial mass Mv, radius and 
velocity K are thus 

GMv 



(4) 



i?v ' (47r/3)i?3 
which lead to the approximate virial relations 

V200 - mU\i + z)'/\ Rwo-Mli\l + z)^\ (5) 

where M12 = M^/W^Mq, F200 = K/SOOkms"!, 
i?ioo = i?v/100kpc, and (1 + z)z = (1 + z)/3. As can 
be verified based on eq. (jA6p in eq. ([5|) is accurate 
to a few percent at z > 1. 

Based on eq. ^ and eq. jT]) , the halo crossing time 
tv scales with the cosmological time t. 



(6) 



M 
M 



■.M^,ii + zr. 



(7) 



In the EdS regime fj, -> 5/2. With ^ = 0.14 and the 
appropriate value of the normalization factor s, this ap- 
proximation for the average of M(z) was found to be a 
good match to cosmolo gical N-body simu lations to bet- 
ter than 5% for z > 1 (jNeistein fc Dekelll2008bl) . while 
it becomes an underestimate of ~ 20% at z = 0. Sim- 
ilar fitting formulae with slight variations in the values 
of the parameters, valid in differe nt ranges of mass and 
redshift, were proposed by others ([Fakhouri fc Mal[2008l : 
iGenel et al.ll2Q08ll2010D . 

The power of p = 5/2 can be simply understood 
from the following scaling argument based on the Press- 
Schechter (PS) and Extended-PS (EPS) approxima- 
tions o|_graAatati£M cosmol- 
ogy (|Press fc Schechteilll97'ilBond et al.lll99l[ ). A key 
element in the PS formalism is a self-invariant time 
variable, u! cx D{a)~^, where D{a) is the growth rate 
of linear density perturbations (eq. (jAlOP ). The self- 
invariance means that the growth rate of halo mass with 
respect to uj is independent of lu, namely dM/duj = 
const. ^ which implies M oc In the EdS regime, where 
D{a) cx a and a oc t^/^, this gives M cx a~^/^. We 
therefore use p = 5/2 for our toy model at z > 1. 
At lower redsh ifts, a better fit can be obta ined with 
H 2.4 2.2 (iNeistein. van den Bosch fc De kel 2001 
iNeistein fc Dekelll2008bf ). 

The small power f3 reflects the log-slope of the fluc- 
tuation power spectrum. The value 13 ~ 0.14 fits well 
the Millennium simu lation merger trees for M y in the 
range 10"-10"Mo (jNeistein fc Dekel2008b[ ) . Accord- 
ing to EPS, it should be /3 = (n -I- 3)/6, where n ~ —2 
is the power index, P{k) oc fc", on the corresponding 
scales. This implies that /3 should be even smaller for 
smaller haloes. Since /? is rather small, we approximate 
/3 = in the toy model addressed in this paper. 

The normalization factor s is the specific accretion 
rate into a halo of — at z = 0, or the 

inverse of the corresponding accretion timescale, s = 
''"in (Ell)- As described in ijAl the value of s for haloes 
of Mv = lO^^A/0 has been estimated earlier to be s ~ 
0.030 Gyr~^, by a fit to the systematic mass growth 
A'I{a) in the Millennium cosmological simulation, scaled 
to the A CDM cosmological pararn eters assumed here 
f following fNeistein fc E)ekell l2Q08bl) . In particular, the 
quoted value of s is for as = 0.82, and it scales as s cx 



3.2 Toy Model: Accretion 

We use the general term "accretion" to refer to the total 
infiow into a distinct halo or its central galaxy, including 
all the mass in dark matter, gas and stars and not dis- 
tinguishing between smooth and clumpy components, 
i.e., including all mergers. 

As reported in the Appendix, especially ^A41 the 
average specific accretion rate of mass into haloes of 
mass M at z can be approximated by an expression of 
the form 



According to eq. ([7]), with s ~ 0.030 Gyr~^, a halo 
of at z = 2 accretes baryons at a rate Mi,,ac — 

80/bo.i7 M© yr~^, where /bo. 17 is the universal baryon 
fraction in units of 0.17. 

In a simple toy model that we address here, valid 
for massive galaxies at z > 1, we ignore the weak mass 
dependence in the average specific accretion rate, /3 — ^ 
0, and simplify eq. ([7]) to 



s ~ 0.030 Gyr" 



(8) 
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Figure 2. Cosmological accretion of total mass: specific mass 
inflow rate M/M through the virial radius as a function of time 
(expansion factor a and redshift 2). Shown at each time step 
are the median (thick black), the average (magenta), and the 
average of the log within the 90% percentiles (green) over the 
sample of simulated galaxies. The error bars estimate the error 
of the mean, and the shaded area marks the 68% percentiles. 
Shown in comparison is the toy-model prediction, eq. JHJ with 
s = 0.030 Gyr~^ (thick smooth red). The toy model provides a 
reasonable fit over the whole redshift range 2 = 4 — 1 for the 
median and for the averages. 
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Figure 3. Cosmological accretion of total mass: growth of virial 
mass. Shown arc the median (thick black), the average (magenta), 
and the average of the log within the 90% percentiles (green) over 
the simulated galaxies. The mass of each galaxy M(z) has been 
scaled before stacking by median[Af(2: = 2)]/M{z = 2) (see text). 
Shown in comparison is the toy model prediction, eq. ^ with 
a = 0.79, normalized like the simulations at 2 = 2 (thick smooth 
red) . The toy model is a reasonable approximation over the whole 
redshift range. Also shown is the toy model with q = 1.1 (dashed 
red), which fits the median in the range 2 = 4 — 2. 



Figure 4. Cosmological accretion of total mass: mass inflow rate 
M. Shown are the median (thick black), the average (magenta), 
and the average of the log within the 90% percentiles (green) 
over the simulated galaxies. M{z) of each galaxy has been scaled 
before stacking by mcdian[M(2 = 2)]/M(2 = 2) (sec text). Shown 
in comparison is the toy model prediction, eq. I|10|l with a = 
0.79 and s = 0.030 Gyr~^, normalized like the simulations at 
2 = 2 (dashed red). The model reproduces the stacked simulation 
results to 0.1 — 0.2 dex. 



This can be simply integrated to a growth of halo mass 
as a function of z, where the mass at some fiducial red- 
shift zq is given to be Mq, 



My = Mo e 



— a (2 



(3/2) sti ~ 0.79. 



(9) 



This functional form was indeed found to be a good 
fit to the halo growth in earlier cosmological jV-body 
simulations of lower resolution (jWechsler et alj 120021 : 
Neistein k Dekel 20081s). 

Note from eq. ^ and eq. ^ that the accretion rate 
into a given halo as it grows is 



M{z) ~ s Mo e"" (1 + z)5/2 . 



(10) 



This average accretion rate into a given halo as it 
evolves does not vary much in time over an extended 
cosmological period. For a = 0.79, this rate has a max- 
imum at z = 2.5/ a — 1 ~ 2.2, and it varies by less than 
a factor of 2 in the range z ~ 0.3 — 5. The maximum 
average baryon accretion rate in the history of a halo 
of mass Mo = 2 X IO^^Mq today, similar to the Milky 



Way, is Mn 



33/bo.i7^0yr" 



While the expressions derived above are for the to- 
tal accretion rate dominated by dark matter, we sus- 
pect that the same expressions for the specific accre- 
tion rates could be valid for the specific accretion rate 
of baryons into the virial radius. This should be true 
when the baryons follow the total-mass inflow with a 
constant baryonic fraction, and as long as we tenta- 
tively ignore the baryonic mass loss from the haloes. 
Note that strong outflows can in principle make the 
net baryonic accr etion rate smaller than the total ac- 
cretion rate (e.g. iFaucher- Giguere" Keres fc Mai 120111 : 
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Ivan de Voort et all l2011a^ , but the specific accretion 
rate may remain the same. 

3.3 Simulations: Accretion 

3.3.1 Total accretion 

In Figs. [5] to m we address the total mass inflow rate, 
dominated by the dark matter component and includ- 
ing the baryons, gas plus stars. The accretion rate 
M through a spherical boundary of radius R during 
a timestep At is taken to be the difference between the 
masses encompassed by the sphere of radius R at the 
two snapshots defining the beginning and the end of the 
timestep, divided by At. In most cases R is either R^ 
or 0.1i?v, and in each of the snapshots we use the ac- 
tual value of i?v at that time. We verify in ^4.2i Fig.fTTl 
that using a fixed radius docs not make a significant 
difference at z > 1. An alternative way to compute the 
instantaneous M in a given snapshot would have been 
to sum over the cells or particles of mass rrii and ra- 
dial velocity Vr^i within a shell of thickness D about R, 
namely M = rrii Vr.i/ D. As a sanity check, we have 
verified that when averaged over a timestep this gives 
similar results to the former method. 

We highlight in each figure the median over the 
galaxies within timesteps of Aa = 0.02, and also show 
the linear average, as well as the average of the log 
within the 90 percentiles. Because of the lognormal na- 
ture of the distribution, the average of the log tends to 
be similar to the median, while the linear average tends 
to be larger. These stacked quantities are compared to 
the corresponding toy-model prediction. 

Figure [2] shows the specific total mass infiow rate 
M /M through the virial radius. The simulation results 
are compared to the toy- model prediction, eq. (jH]), with 
s = 0.030 Gyr"\ There is no other free parameter in 
the toy model for M jM. This comparison directly tests 
the predicted asymptotic systematic time dependence, 
[i — b/2, with the approximation that the specific infiow 
rate is independent of mass within the relatively narrow 
mass range spanned by the current sample. This mass 
independence allows a straightforward stacking of the 
galaxies of different masses without any scaling. This 
makes the specific accretion rate the most robust quan- 
tity for the comparison of the model with the simula- 
tions, and for determining the best-fit normalization, 
the parameter s. We note, however, that a success in 
predicting the systematic evolution oi M/M does not 
guarantee an accurate match to the evolution of M and 
of M, which could be both off by the same multiplica- 
tive factor. Another mismatch may occur because the 
average of M /M is not necessarily the ratio of the av- 
erages of AI and M. 

One can see in Fig. [2] that, in general, the model 
provides a good fit to the systematic rcdshift depen- 
dence in the simulations, and that s — 0.030 Gyr^^ 
gives a proper normalization. The model matches the 
median (and log average) of the simulations well in the 
redshift range 2 = 4 — 2, where our sample is complete 



and where we focus our analysis. The model is a slight 
overestimate of the median in the range z = 1.8 — 1.3, 
with deviations at the 2(7 level. As expected, the linear 
average is somewhat higher than the median at every 
time, and it is approximately matched by the model 
throughout the whole redshift range z = 4 — 1, though 
the average is slightly above the model at z = 2.6 — 2.0. 

Figure [3] shows the corresponding growth of total 
virial mass in time. In order to stack all the galaxies, 
the mass of each galaxy at z has been scaled by the in- 
verse of the ratio of its mass at z = 2 to the median mass 
at z = 2, namely M{z) is multiplied by mcdian[M(z = 
2)]/M(z = 2). For the stacking at a given z, the aver- 
age (M(z = 2)) is computed for the sample of galaxies 
that is averaged over at z, namely the whole sample 
of 27 galaxies at z ^ 2, and the gradually diminish- 
ing subsample from z = 2 to 1. The simulation results 
are compared to the toy model prediction, eq. © , with 
a = 0.79, corresponding to s = 0.030 Gyr~^. The model 
is normalized at zq = 2 to match the simulations' me- 
dian, Mq = median[M(z = 2)], as defined above. This 
way of stacking yields by construction no scatter and 
perfect fit of the model at z = 2; the comparison thus 
tests the success of the model at redshifts away from 
z = 2. We learn that this model is a reasonable ap- 
proximation for the total mass growth over the redshift 
range studied. The model overestimates the median 
by ^ 0.2 dex at z = 4 — 3, but it is less than 0.1 dex 
above the linear average. The figure demonstrates that 
a higher value of a (and s) can fit better the average or 
the median in the range z = 4 — 2, where the simulation 
results are more reliable. 

Figure [4] shows the corresponding absolute total 
mass infiow rate M at R^. Motivated by eq. ([TU)) where 
AI (X M, the value of M for each galaxy at every z 
has also been scaled before stacking by median[M(z = 
2)]/M(z = 2) (similar to the scaling applied to M(z)). 
The simulation results are compared to the toy model 
prediction, eq. ([TU]), with s = 0.03 Gyr~^ and a = 0.79. 
The toy model is again normalized at zq = 2 to Mq ~ 
median[M(z = 2)]. Indeed, the value of M docs not 
vary much in time during the evolution of a halo within 
this redshift range, because the growth of halo mass 
compensates for the cosmological decline of specific in- 
fiow rate at a given mass. The gradual diminishing of 
the sample between z = 2 and z = 1 is reflected in the 
shape of the model curve in that region. The model re- 
produces the general trend with time seen in the stacked 
simulations. It overestimates the median by ~ 0.1 — 0.2 
dex at z = 4—2.7, and even by more at z '--^ 1.4—1.3, but 
it matches well the linear average over the whole range 
(except at z < 1.3 where the sample is very small). 

Our results for the total accretion at R^ confirm 
that the toy model in its simplest form provides a very 
useful approximation for the characteristic values of the 
three quantities, M/M, M and M, over the whole red- 
shift range. We see that the simplest toy model enjoys 
different levels of success for the different quantities at 
different redshifts. Somewhat different values of s and 
a may provide the best fit to the median and to the 
average, to each of the three quantities, and at differ- 
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(Mbar/Mbar)/model log [(Mbar/Mi,ar)/moclel] log [(Mba,/Mb3,)/model] 

Figure 5. Baryon accretion: Distribution of specific accretion rate of baryons, M/M, relative to the model prediction, as measured 
at iJv (left linear, middle log) and at O.liJv (right). All snapshots of all galaxies are included, averaged within time steps of Aa = 0.02. 
The distributions are close to lognormal, as can be seen from the lognormal functional fits, with the associated logio mean {/iio) and 
standard deviation (crio) quoted. The dashed Gaussian in the right panel is the fit at Rv from the middle panel. The scatter represents 
variations among the galaxies and along the history of each individual galaxy (see Fig.O. Marked by vertical bars are the median and 
linear average. The average of the log within the 90% percentiles almost coincides with the median. The gray and yellow shaded areas 
denote the 68 and 90 percentiles about the median. The negligible tail of negative values (excluded from the log plot) corresponds to 
two rare cases. The few percent tail of snapshots with high accretion rates corresponds to mergers. 
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Figure 6. Baryon accretion: specific inflow rate and penetration. 
Shown are the median (thick solid) and the average (thin dashed) 
of M/M through Rv (black) and through O.lRv (green). At Rv 
it is analogous to the total rate in Fig. [2] The model prediction 
for i?v, eq. ((Sjl with s = 0.030 Gyr~^ (smooth solid red), is in- 
deed a good approximation at Rv, but it also provides a good 
approximation for the specific accretion rate at O.l/Jv 



Figure 7. Baryon accretion: specific infiow rate at Rv for four 
individual galaxies (that enter the average in Fig. Ell. The black 
curve second from bottom is normalized properly, and the other 
curves are shifted by 2 dex relative to each other. The simu- 
lations are compared to the toy model prediction, eq. ^ with 
s = 0.030 Gyr~^ (dotted). This figure illustrates the scatter due 
to variations among the galaxies and due to variations along the 
history of each individual galaxy. 



ent redshift ranges. There are several reasons for this. 
First, the lognormal nature of the distribution of M 
makes the average larger than the median. Second, the 
three noisy quantities are not linearly related to each 
other, and therefore their averages do not simply re- 
late to each other. Third, the toy model in eq. ([8]) and 
eq. ([9]) slightly deviates from the more accurate approx- 
imation eq. (O, where /3 ~ 0.14, which itself is only an 
approximation. 



3.3.2 Accretion of baryons 

The baryonic inflow through the virial radius in the sim- 
ulations consists of stars and gas in comparable frac- 
tions, with the stellar fraction increasing with time. 
Most of the stars are in small merging galaxies, associ- 
ated with mini-minor m ergers of mass ratio smaller than 
1:10 (|Dekel et al.ll2009D . These mergers do not have a 
substantial global dynamical effect on the disc — they 
tend to join the disc as part of the cold streams and 
grow the cold disc component that develops the vio- 
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Figure 8. Baryon accretion: growth of mass within i?v, analo- 
gous to the total mass in Fig. [3] The toy model with the adopted 
value of a eq. ||9ll is a reasonable approximation over the whole 
redshift range. Also shown is the toy model with a = 1.1 (dashed 
red), which fits the median in the range 2 = 4 — 2. 



lent disc instability. On the other hand, the incoming 
gas component provides much of the fuel for new star 
formation. Recall that the stellar fraction at z ~ 2 is 
considered to be an overestimate in the current simu- 
lations, by a factor of order two, because of the high 
efficiency assumed for star formatio n at high redshift 
(jCeverino. Dekel k, BournaudI l2010l) . We therefore fo- 
cus on the more robust total baryonic accretion, and do 
not seek an accuracy of better than a factor of two when 
dealing with the gas accretion alone, e.g. in Sj6]and ^ 
Figures [5] to [8] address the accretion of baryons 
through the virial radius. The calculations for the bary- 
onic accretion are analogous to those described in i )3.3.1l 
for the accretion of total mass. The first three figures 
refer to the specific rate, where Fig. [5] shows the distri- 
bution over all snapshots. Fig. [6] shows the systematic 
time evolution over the simulations compared to the toy 
model prediction eq. dS]), and Fig. [7] shows the same for 
four individual galaxies. Then Fig. [prefers to the baryon 
mass growth in the simulations compared to the model 
prediction eq. 

A comparison of Figs.|n]and[H]to Figs. [2] and |3] shows 
a strong similarity, implying that the baryons follow the 
dark matter in its streaming into the virial radii. This 
implies that the outflows through the virial radii are 
indeed rather small, at least on average in this suite of 
simulations. 

In Fig. [5] we learn that the distribution of specific 
baryonic accretion rates at i?v resembles a lognormal 
distribution. This is very similar to the distribution of 
the total accretion rates, which has been seen earlier 
to resemble a lognormal distribution ([Neistein fc Dekel 
l2008a[) . The median and averages of the ratio oi M /M 
to the model prediction are close to unity, indicating 
that the toy model provides a good approximation for 
the systematic behaviour. The median and the very sim- 



ilar average of the log are at 0.92, and the linear average 
is 1.23, somewhat larger as expected from a lognormal 
distribution. This is an argument for using the median 
or the average of the log to represent the systematic 
behaviour of the stacked sample. The scatter of ±0.33 
dex indicates that the specific accretion rate in different 
galaxies or at different times can deviate significantly 
from the median or averages. Small negative M values, 
namely net outflows, are measured in two rare cases, 
reflecting a combination of satellites moving out of the 
halo and gas outflows. The extended tail at large accre- 
tion rates corresponds to mergers. 

Figure [5] shows that, similar to the case of total 
accretion, the toy model of eq. ^ provides a good 
flt to the stacked simulation results at R^. With s = 
0.03Gyr-\ it is an excellent match to the median at 
z = 4 — 2. The model is a slight overestimate of the 
median at z = 1.8 — 1.0, at the level of 0.1 — 0.2 dex, 
or 2-sigma in terms of the error of the mean in a few 
points. With this normalization, the average is slightly 
higher than the model at 2 = 4 — 2, and is matched well 
at z = 1.8 - 1.2. 

Figure [7] shows the speciflc baryon accretion rate 
at i?v for four individual galaxies. It shows separately 
the scatter due to variations among the galaxies and 
the scatter due to variations along the history of each 
individual galaxy. Although some galaxies seem to de- 
viate systematically from the model throughout their 
histories (e.g. the second from bottom case tending to 
be lower than the model at most times), most galax- 
ies fluctuate about the model prediction. The typical 
scatter along the histories of individual galaxies is com- 
parable to the variations from galaxy to galaxy. 

Figure [5] addresses the baryon mass growth within 
i?v In analogy to Fig. [31 the galaxies are scaled before 
stacking to the same value of the median at z = 2. The 
model overestimates the median by 0.2 dex at z ^ 4, 
but it overestimates the average only by 0.1 dex at that 
redshift. A higher value of a (and s) matches the growth 
curve better in the range z = 4 — 2. 



4 PENETRATION TO THE INNER HALO 
4.1 Toy model: Penetration 

We expect, based on the reasons outlined below, that at 
z > 1 a very large fraction of the baryons that enter the 
halo at the virial radius efficiently penetrate into the 
galaxy at the halo centre. Below a critical halo mass 
of ^ 1O^^M0, a stable virial shock encompassing an 
extended hot gas medium cannot be supported because 
the radiative cooling rate is too high (iBirnboim fc Dekell 
l2003tlKeres et al.ll2005tlDekel fc I3irnboimll2006f ) . so the 
instreaming of cold gas is unperturbed. A stable virial 
shock is likely to develop in more massive haloes, and 
it may shut down the cold accretion in sufficiently mas- 
sive haloes at low redshifts, but at z > 1, most of the 
accreted gas is expected to be in cold streams that pen- 
etrate thr ough the hot halo me dium deep into t he in- 
ner halo (lOcvirk Pichon fc Tev ssier 2008; Deke l et all 
[2009HKerei et al.ll2009D ." as illustrated in Fig. [5] (based 
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Redshift z 

Figure 9. Predicted penetration of cold gas streams into 
tlie lialo centre as a function of halo mass and redshift . This 
schematic diagra m (reproduced from iDekel fc BirnboimI |2006| : 
iDekel e t al.||20 0"9|). is based on analytic spheri-sym metric calcula- 
tions llDekel fc Birnboimll2006l: fPekel et al.ll200g|) . While a virial 
shock is expected to be present at and above Mv ~ 10^^ Mq, at 
high redshift the gas is expected to penetrate through it in cold 
streams along the filaments of the cosmic web. 
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Figure 10. Baryon accretion: inflow rate and penetration. Shown 
are the median (thick solid) and the average (thin dashed) of 
M through iJv (black) and through O.liJv (green). At Rv it is 
analogous to the total rate in Fig.|4] The model prediction for Rv, 
eq. ([Toll with a = 0.79, s = 0.030 Gyr~\ and normalized like the 
simulations at z = 2 (smooth solid red), is a good approximation 
at Rv . The penetration to the inner galaxy is ~ 50% at 2 = 4 — 2, 
and higher at z = 2 — 1. 



on'Dck cl fc BirnboimI |2006( ). At high redshift, in mas- 
sive galaxies that represent high-sigma peaks in the den- 
sity fluctuation field, the gas density in these streams 
is high compared to the mean gas density in the halo, 
as they follow the narrow dark-matter filaments of the 
cosmic web. The high density enhances the radiative 
cooling rate in the streams and prevents the devel- 
opment of pressure that could support a stable virial 
shock in the streams. As a result, the streams pene- 
trate through the halo with a mass inflow rate that is 
rather constant with radius (IDekel et al.ll2009l ). They 
flow in at a constant inflow velocity that is compara- 
ble to the virial velocity, and they keep cold at lO^K 
while dissipating the gained gr avitational energy into 
cooling Lyman-alpha radiation (iDiikstra fc Loebll2009l : 
iGoerdt et aLll2010l: iFumagalh et al.ll201lN H 



4.2 Simulations: Penetration 

Figure |6l which has been discussed above concerning 
the specific accretion rate of baryons at , also shows 
the specific accretion rate of baryons into the inner halo 
at 0.1i?v The stacked simulation results are compared 
to the model prediction for the accretion rate at i?v at 
that time, eq. (|S|). The right panel of Fig. |S] shows 
the distribution of specific accretion rate of baryons at 
0.1i?v, over all the snapshots of all the galaxies, normal- 
ized to the model prediction. The distribution at 0.1i?v 



^ Note that SPH simulations, which tend to underestimate the 
dissipat ion rate in the cold streams, may yield somewhat different 
results jFaucher-Giguere et al] |2010| : Ivan de Voort et alj l2011bl : 
iNelson et alJl2013lV 
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Figure 11. Baryon accretion: inflow rate and penetration at 
fixed radii. Same as Fig. 1101 except that the evolving i?v and 
O.liJv are replaced by fixed radii 100 kpc and 10 kpc, respectively. 
There is no qualitative change from Fig. 1101 at 2 > 2, indicating 
that the artificial "accretion" due to the growth of R^ is not a 
major part of the actual inflow rate into the central galaxy at 
high redshift. 



should be compared to that shown in the left panels for 
Rv, and to the Gaussian fit that characterizes it. We 
learn that the specific baryon accretion rate at i?v and 
at 0.1 i?v arc rather similar, and well approximated by 
the toy model. In fact, the model is a good match to 
the median of the simulations at 0.1 i?v over the whole 
redshift range. 
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The similarity between the specific accretion rate 
at 0.1i?v and at does not necessarily imply 100% 
penetration through the halo, as both M and M for the 
baryons could decrease between i?v and 0.1 i?v in a simi- 
lar way. We conclude that eq. ([5]), with s = 0.030 Gyr"^, 
is also a useful approximation for the specific accretion 
rate of baryons into the galaxy itself. This is for the 
distribution among galaxies of halo masses that range 
from below 10"Mq to above IO^^Mq in the redshift 
range 2 = 4—1. 

Figure [10] shows the absolute inflow rate M for the 
baryons at and at 0.1i?v) and the corresponding toy 
model prediction at i?v^ eq. (fTO)) . The scaling to the 
average at z = 2 before stacking the simulated galaxies 
and the corresponding normalization of the toy model is 
the same as described in ^13.3.11 Wc first note in Fig. [10] 
that the simplified toy model, eq. (|TU|) . is a reasonable 
fit to the inflow rate of the simulations at i?v in the 
range z = 4 — 1. The model overestimates the median 
by 0.2 dex or less, and it fits the average to within 0.1 
dex (except at z < 1.3 where the sample is very small). 
Indeed, the value of M does not vary much in time dur- 
ing the evolution of a galaxy within this redshift range, 
because the mass growth compensates for the decline 
of specific inflow rate, as for Fig. [4] This fit is comple- 
mentary to the good fits provided by the model to the 
specific inflow rate and to the mass growth. Figs. [5] and 

m 

Figure [10] provides a new insight into the penetra- 
tion of baryons through the halo into the inner galaxy. 
One reads that, in the range z =4 — 2, the baryonic 
inflow rate at 0.1 i?v is typically ~ 50% of the baryonic 
inflow rate at R^. In the range z = 2 — 1 the penetra- 
tion seems to be more efficient, though the uncertainty 
is larger because the sample is smaller. This massive 
penetration is in general agreement with earlier esti- 
mates base d on the MareNostrum simulation with l ower 
resolution ijPekel et al.ll2009l: iDanovich eral]|2012[ ) and 
with the high SFR observ ed in massive galaxies at z ^ 2 
(iGenzel et al.ll2008L l201lD . If the SFR follows the gas ac- 
cretion rate into the disc (see ^ , the SFR can in prin- 
ciple be a fixed significant fraction of the instantaneous 
overall baryon accretion rate into the virial radius. Of 
course, internal processes within the disc may suppress 
the SFR at early times (z > 3, sa y) and less massive 
haloe s (Mv < IO^^Mq, say) (e.g. iKrumholz fc Dekell 
|2012|) . and thus accumulate gas for higher SFR in the 
more massive galaxies at later times (e.g., z ~ 2). 

One should recall the caveat associated with the sig- 
nificant fraction of stellar mass in the inflow through the 
halo into the central galaxy. This may be partly respon- 
sible for the more efficient penetration seen at z < 1.8, 
and it reminds us not to assume accuracy of better than 
a factor of two when addressing the incoming gas and 
the resultant star formation rate. 

The radii of spheres through which the inflow rate 
has been considered so far were the virial radius for 
the outer halo, and a fixed fraction of it (0.1i?v) for 
the galaxy at its centre. This is in order to allow us 
to consider self-similar radii that follow the halo and 
the galaxy as they grow. However, the growth of R^ 



in time (which can be deduced from the virial rela- 
tions, eq. ([S|), and the toy model for halo mass growth, 
eq. ((HI)) is responsible for part of the mass growth within 
i?v, which may not be associat ed with actual inflow 
()Diemer. More fc Kravtsov]|2013[) . In order to quantify 
this effect, Fig.fTTjshows the inflow rate through spheres 
of fixed radii, 100 kpc and 10 kpc, replacing i?v and 
0.1i?v of Fig. [TU] A comparison between the two fig- 
ures indicates no significant differences between the ab- 
solute values of M at z = 4 — 2, and similar conclusions 
regarding the penetration efficiency. Similar results are 
obtained for the specific accretion rate. We conclude 
that, at least at z ^ 2, the vast majority of the inflow 
measured through R^ or a fixed fraction of it is true in- 
flow, and therefore proceed with the analysis using the 
radii that grow self-similarly in time, i?v and 0.1i?v 



5 DISC SIZE 

5.1 Toy Model: Disc Size versus Virial Radius 

We assume that the characteristic disc radius scales 
with the halo virial radius via a contraction factor A, 

i?d = Ai?v. (11) 

We also assume that the characteristic circular velocity 
of the disc scales with the virial velocity, 

Va = vV^. (12) 

Given that the halo is roughly an isothermal sphere, 
the value of v is expected to be of order unity, and 
we sometimes adopt v ~ a/2, which is a good ap- 
proximation for the Milky Way. The value of A can 
be estimated in the common model where the galactic 
disc is assumed to form by di ssipative gas contraction 
within the dark-inatter halo (iFall fc Efstathioul Il980l : 
iMo. Mao fc Whit3[l998l: iBullock et al.ll200lD . If one as- 
sumes that the original specific angular momentum of 
the gas, J, is similar to that of the dark matter in the 
virializcd halo, and if j is conserved during the gas con- 
traction, then 

\~y/2v-^\', (13) 

where A' = j / (^/2R^V^) is the halo spin parameter as 
defined bv IBullock et al.l ()2001l) . Based on N-body sim- 
ulations and tidal-torque theory, the halo spin param- 
eter is assumed to have a constant avera ge value inde- 
pendent of mass or time, A' ^ 0.035 fe.g. IBullock et al] 

l200lt) , implying that A is constant and of a similar value. 
The disc dynamical time is then 

id = ^ ^ v~^\ ^ ~ 0.0071 v-%.05 t , (14) 

where A0.05 = A/0.05, and where we have used the rela- 
tion of the virial crossing time to the cosmological time 
t, eq. (O. 

Note that the distribution of spin parameter among 
haloes is lognormal with a standard deviation of half 
a dex (iBuUock et al.|[200ll ). This would translate to a 
large scatter about the average i?d of eq. (|TT|) . 
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Figure 12. Ratio of disc radius to virial radius, A. Shown is the 
median over the simulations (black curve) and the 68 percentiles 
(shaded area) for the snapshots where there is net inflow, assumed 
to undergo VDI. The dotted line in the range z = 4 — 3 reflects 
the large uncertainty in the disc radius in that regime. The ratio 
is in the ballpark of, and slightly higher than, the expected value 
for the halo spin parameter, validating eq. The average is 

slowly declining in time. 



Using eq. ((TT|) and eq. ([5]), one can derive a useful 
expression for the baryonic surface density in the disc, 



~ 1.3xlO^M0kpc-2md,o.iA^.^5Af^/'(l + z)2 , (15) 

where Md is the cold mass in the disc, and md,o.i = 
md/0.1 refers to the baryonic mass in the disc relative 
to Mv. 



fer to de-emphasize the results that explicitly depend 
on the disc radius at z > 3, e.g., in Fig. [TSland Fig. [T71 
below. 

Figure [T^ shows the average and 68% scatter of the 
ratio i?d/i?v over the simulated galaxies. In eq. ([TT]) . a 
value of A ~ 0.05 crudely fits the 1-sigma range of val- 
ues from the simulations over the whole redshift range. 
However, there seems to be a systematic trend with red- 
shift, from a median of ~ 0.06 at z = 3 — 2, through 
0.05 at z = 2 - L5 down to 0.04 at z = 1. The in- 
dicated values of A are somewhat larger than what is 
implied from the spin parameter of dark matter haloes 
as estimated from cosmological N-body simulations, es- 
pecially if V is larger than unity. T his is consistent wit h 
observational indications at z ~ 2 (jGenzel et al.|[2006h . 
Indeed, the specific angular momentum of the baryons 
in the high-redshift discs is likely to be higher than that 
of the dark-matter haloes (a) because of the way angular 
momentum is transpor ted into the disc by cold streams 
from the cosmic web (iKimm et al.l 201 It Pichon et al.l 
120111 : iDanovich etHI l2012t IStewart et al.l l2013| ). and 
(b) because outflows tend to pre ferentially carry away 
low- angular-mom ent um material ([Mailer fc Dekelir2002l : 
iBrook et al.ll201l( l. We take the toy model for the disc 
radius and the associated estimate of dynamical time, 
eq. (fTTj) and eq. (fT4)) . to be accurate to within a factor 
of two. 

Using abundance matching to determine halo virial 
radii, iKravtsovl (|2013l ) finds that the half-mass radii of 
today's galaxies of all morphological types, i?haif, are 
related to their today's halo radius i?2oo as i?hait — 
0.015i?2oo- If the galaxy has formed roughly in its 
present size at z ~ 2, when the halo was a few times 
smaller in radius, the implied ratio then is consistent 
with our finding for A in Fig. 1121 



5.2 Simulations: Disc Size versus Virial Radius 

The disc ra dius Rd is determined in the simulations as 
described in iMandelker et ahl (|2013| ). In short, the disc 
is modeled as a cylinder whose z axis coincides with 
the spin axis of the cold gas within it (T < 1.5 x 10'*K, 
which is typically 97% of the gas) . We start with a large 
cylinder of radius 0.15i?v and half-height Ikpc, and it- 
eratively converge on the cylinder of radius Rd that con- 
tains 85% of the gas mass within the large cylinder. This 
determination of the disc radius is commonly in good 
qualitative agreement with the radius one would have 
estimated by visual inspection of the face-on surface 
density of the cold gas. Note, however, that this is only 
a crude approximation for the radius of the whole disc 
including its stellar component, which could be domi- 
nant. 

The disc radii as determined from the simulations 
at z > 3 are rather uncertain and sometimes ill-defined. 
This is because the discs in these early times tend to be 
small and highly perturbed as the timescale for mor- 
phologically damaging mergers is comparable to and 
shorter than the disc dynamical time. We therefore pre- 



6 MASS CONSERVATION: STEADY-STATE 
OF SFR AND GAS MASS 

6.1 Toy Model: Steady State 

Assuming mass conservation, the gas mass in the galaxy 
varies subject to a so urce term and a sin k term (analo- 
gous to a "bathtub" . iBouche et al.|[2010t ). e.g.. 



A4 = M,, 



Ms, 



(16) 



where Mg^ac is the accretion rate of gas, and Mgf is the 
star formation rate. We temporarily ignore potential gas 
outflow from the galaxy (to be incorporated later). We 
assume that the SFR density obeys a universal local vol- 
umetric law, psf = CsfPg/iff with Egf 0.02 a constant 
efficiency and pg th e density of molecular gas in th e 
star-forming region l|Krumholz. Dekel fc McKe3l2012[ ). 
Then, under certain circumstances, the overall SFR in 
the galaxy can be crudely assumed to be proportional 
to the total gas mass, 

Msf = , 

where t^i = tff/est- 



(17) 
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If Mg^ac and t^j^ vary on a timescale longer than 
Tgf , the simple solution of eq. (|16p with eq. (|17p is 

Mg = A/g,ace-*/"='' , Afg = A'/g,acrsf(l-e^*/"='')- (18) 

After a transition period of order Tgf, the solution re- 
laxes to a steady state solution 



, M„ ~ Mg,acTsf 



(19) 



The SFR sink term adjusts itself to match the external 
source term Afg.ac- 

The value of Tgt can be related to the global disc 
crossing time by iff = /sfid, with f^i ^ 0.5, assum- 
ing that star formation occurs in regions (partly bound 
clumps) where the overdensity is a few with respect to 
the mean density in the disc. Using the estimate for 
from eq. ([13]), we obtain, for fsf ^ 0.5 and egf ~ 0.02, 



r,f = e-i/,ftd-0.17t. 



(20) 



At z = 2, where the Hubble time is i ~ 3.25 Gyr, we 
have Tsf - 500 - 600 Gyr. 



The validity of the solution eq. ([18)) to eq. ([161) de- 
pends on the timescales for variation of Mg.ac and of r^^ ^ 
compared to Tsf . Using eq. ([5]) and eq. ([H wc obtain 



(A^g.ac) 



MMg,acM| |(2/3)aa-i-5/3| 



(21) 



With a = 0.79 in eq. ([9]), we get tvar(A^g,ac)/t ^ 
1.0, 2.3, 12, 1.6 for z = 4, 3, 2, 1 respectively. This means 
that tvar(A/g,ac) ^ ^ throughout the redshift range 
z = 4 - 1, so, based on eq. ([^Dj), ivar(-/W^g,ac) > ^sf, 
as required. Similarly, using eq. pop. the timescale for 
variation of is 



^var ('Tjf ) 



'sf 



(22) 



namely |ivar(Tgf^)| S> Tgf, as required. We thus expect 
the steady-state solution to be valid through the whole 
period of interest here, and to be approximately valid 
out to z 10 and earlier. 

A continuity equation of a similar nature, 
and its steady-state solution, were used with a 
variety of sink terms in several recent studies 



(iDek cl. Sari fc Ceverinol 
Krumholz & Dckej l20lf 



2009t_ 



20101: 



2012t iGenel. Dekel 



.Bouche et al.l . _ 

ICacciato. Dekel fc Genell 

Cacciaton20"l2l: lFeldmannll20100 . 



The sink term in eq. (|16p can be generalized to include 
other sink terms. For example, if there is an outflow at 
a rate A'/out = /outA^sf , the second term of eq. ([TO]) is 
simply multiplied by the factor (1 + /out)- Similarly, 
if eq. (jl6p is applied only to the disc component, 
instability-driven mass inflow within the disc adds 
another sink term of a similar form, see S[7| In these 
cases, Tsf in eq. ([T5)) is replaced by a smaller timescale 
r, which makes the steady-state solution an even better 
approximation. 
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Figure 13. Overall steady state: Shown are SFR and Mgas 
within the sphere of radius 0.1-Rv, as well as the gas accretion rate 
Mgas,acc through that radius, which is the sum of the former two. 
Shown are the medians (solid thick) and the average (dashed thin) 
over the simulated galaxies. The 68 percentiles (shaded area) and 
the errors of the mean (error bars) are shown for A/gas and SFR. 
The simulations evolve about the steady-state solution, eq. I|19|l . 



where SFR= M^, 



and Ms-. 



■ 0, with a slow increase in the 



2:as mass at the level of ~ 10% of the SFR and the accretion rate. 



6.2 Simulations: Steady State 

FigurcsfinitolTSjverifv the validity of the steady state so- 
lution, eq. (|19p . in the simulations. Figure [T3l shows the 
SFR and the rate of change of gas mass, Mgas, within 
the inner sphere of radius 0.1i?v As in ^3.3. 1[ these 
quantities were scaled before stacking to match the me- 
dian mass at z = 2. The gas accretion rate through that 
radius. A/gas, acc, is the sum of the two. We see that the 
median and average of Afgas is indeed small, at the level 
of ~ 10% of the SFR or the accretion rate. This corre- 
sponds to the SFR being close to the accretion rate, 
typically only ~ 10% smaller. 



Figure [T4| shows the ratio of SFR and M^ 



gas, acc 



at 



0.1 i?v for four individual galaxies, compared to the 
steady-state prediction of unity, illustrating the vari- 
ation along the history of each galaxy and among the 
different galaxies. As in Fig. [71 we see that the typical 
galaxy does not tend to show a systematic deviation 
of the SFR from the gas accretion rate into the disc 
along its history. The scatter along the history of each 
individual galaxy is similar in the different galaxies and 
comparable to the variations between different galaxies. 

Figure [TS] shows the distribution of SFR/Afgas,acc 
at 0.1i?v over all the snapshots. It is close to lognormal, 
with a logio mean fXiQ = —0.04 and standard deviation 
(Tio = 0.12. The linear average is only 2% off unity. We 
conclude that the SFR closely follows the gas accretion 
rate into the 0.1 i?v sphere, as predicted by the steady- 
state solution, eq. ([T^ . While the scatter in this ratio 
is only at the level of ^ 30%, the much larger scatter in 
SFR seen in Fig. [T3] indicates large variations in both 
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Figure 14. Overall steady state: SFR versus gas accretion rate 
at O.liiv for four individual galaxies. The black curve second 
from bottom is normalized properly, and the other curves are 
shifted by 2 dex relative to each other. The simulations are com- 
pared to the toy-model steady-state solution, cq. I I19I I. where 



SFR= Mg, 



(dotted). This figure illustrates the scatter due 



to variations among the galaxies and due to variations along the 
history of each individual galaxy. 
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Figure 15. Overall steady state: distribution of the ratio of SFR 
and gas accretion rate at O.liJv AH snapshots of all galaxies are 
included, averaged within time steps of Aa = 0.02. The distri- 
bution is close to lognormal, as can be seen from the lognormal 
functional fit, with the associated logio mean (/iio) and standard 
deviation (crio) quoted. The scatter represents variations among 
the galaxies and along the history of each individual galaxy (see 
Fig. I14|l . Marked by vertical bars are the median and linear av- 
erage. The average of the log within the 90% percentiles almost 
coincides with the median. The shaded areas denote the 68 and 90 
percentiles about the median. The SFR tends to be comparable 
to the gas accretion rate into the central galaxy, with a median 
SFR/Mgasacc ~ 0.91 and an average of 1.02, compared to the 
model prediction of unity. 



accretion rate and SFR between galaxies and/or along 
the history of each galaxy. 



7 INFLOW WITHIN THE DISC 
7.1 Toy Model: Disk Inflow 

Galactic discs at high redshift are expected to de- 
velop a gravitational disc instability with a Toomre 
instability parameter Q smaller than unity (iToomrd 
19641). The Toomre paramete r can be expressed as (e.g. 



Dekel. Sari fc Ceverinol [20091) 



Q 



ttGS 



V2 a 

TV 



(23) 



The radial velocity dispersion a and the surface density 
E refer to the cold component of the disc that par- 
ticipates in the disc instability (see the practical def- 
inition in the simulations in ^17.21 below). We may re- 
fer to it as "gas", but it also includes the cold young 
stars, and one should be careful to specify whether the 
latter are included (see a mu lticomponent analysis in 
ICacciato. Dekel fc Genelll2012|) . The angular velocity ft 
and circular velocity V refer to the characteristic disc 
radius Rd, V = flRd- The constant -\/2 refers to a 
flat rotation curve (and it stands for ■\/2(l -|- 6) where 
V (X r^). The mass fraction 5 is the ratio of mass in 
the cold component to the total mass within Rd, 
Mtot{Rd), the latter including gas, stars and dark mat- 
ter, 

Md 



Mot(i?d) 



(24) 



The instability in high-redshift galaxies is driven 
by the high cold surface density that reflects the high 
gas accretion rate and the high mean cosmological 
density at earlier times. We term this phase Violent 
Disk Instability (VDI), as the associated dynamical 
processes occur on a galactic orbital timescales, as 
opposed to the "secular" processes associated with disc 
instabilities at low redshifts. The unstable disc tends to 
self-regulate itself in marginal instability with Q ~ 1, 
where a represents supersonic turbulence that provides 
the pressure while thermal pressure is negligible. The 
turbulence tends to decay on a timescale comparable to 
the disc dynamical time ([Mac Lowlll999l ). so it should 
be continuously powered by an energy source that could 
stir up turbulence and maintain a at the level required 
for Q ~ 1. In the perturbed disc, which consists of ex- 
tended transient feature s and massive compact clumps 
(jMandelker et al.l |2013[ ). gravitational torques drive 
angular momentum out and cause mass inflow towards 
the c entr e, partly in terms of clump migration (iNqguchil 



19991 Bournaud . Elmegrecn fc Elmegree3 120071: 



Dekel. Sari &: Ceyerino 2009), and partly in terms of 
off-clump inflow (lGammiell200]l : rDekel. Sari fc Ceverinol 
l2009t iBournaud et al.l 1201 it) .This inflow down the po- 
tential gradient from the disc outskirts to its centre 
in turn provides the requ ir ed energy for Q :^ 1 
(iKrumholz fc BurkertI l2010l: ICacciato. Dekel fc Gene! 



I2012t iForbes. Krumholz fc Burkertll2012t iForbes et al 
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|2013[) . The gas inflow rate Mg,i„ can be estimated by 
equating this energy gain and the dissipative losses of 
the turbulence, 



ltd 



(25) 



Here 7 = (2/3)7g ^7$7dis, where 7g is the fraction of 
gas in the inflowing mass, 7$^^ is the energy gain per 
unit mass between the disc radius and the centre, and 
7disid is the turbulence decay timescale. The value of 7 
is of order unity, and it could be as large as a few. We 
thus obtain 



ltd 



1 



27id 



(26) 



assuming Q '--^ 1 in cq. Note that 5 on the right- 

hand side refers to the cold component, including young 
stars, while the left-hand side refers to the gas only. 

Independent estimates based on the mechanics of 
driving the mass inflow by torques yield similar results 
for the inflow rate, to within a factor of two. Exam- 
ples of such calculations arc (a) an estimate of the rate 
of energy exchange by gravit ational clump encounters 



(|Dekel. Sari fc Ceverinori2009l eqs. 7 and 21), and (b) 
an estimate of the angular-momentum exchan ge among 
the t r ansient perturbations in a viscous disc ([Gammiel 
I2OOII: iGenzel et al.l 120081: iDekel. Sari fc Ceverinol 120091 
eq. 24). In Appendix ^A5l we estimate that the disc evac- 
uation rate by clump migration due to the dynamical 
friction exerted by the disc on the clumps is in the same 
ballpark as the above estimates of the inflow rate as long 
as 5 ^ 0.2 — 0.3 as in the steady-sta te solution for VDI 
discs (|Dekel. Sari fc Ceverinoll2009D and in our current 
simulations. Fig. 1181 One should emphasize that the in- 
flow in the disc is a robust feature of the instability, not 
limited to clump migration. Inflow at a rate comparable 
to eq. (p6)) is expected even if bound clumps were to dis- 
rupt by stellar feedbac k in less than a migration time, as 
some t imes assunied ([M urray. Quata ert fc ThompsonI 
l2010t iGenel et al.l l20l3 THopkins ct al.| l2012D . though 
this scenario is not realis tic based on theoretical 
and observational groun ds (jKrumholz fc Dekel I2OIOI : 
IDekel fc Krumholzll2013l) . 

We comment in parentheses that the cosmolog- 
ical streams that penetrate through the halo and 
feed the disc are also potential drivers of turbulence 
(IDekel. Sari fc Ceverinol f200l iKhochfar fc Sili3 I200I 



iGenel. Dekel fc Cacciatol[2012t ) , but this requires strong 
coupling of the streams with the higher-density disc, 
which demands that the streams largely consist of dense 
clumps. Being driven by an external source, this mech- 
anism is not naturally self-regulating. 

Stellar feedback may or may not be a major di- 
rect driver of the turbulence. There is potentially 
enough energy and mo mentum in supernovae, radia- 



tive feedback and winds (iDekel. Sari fc Ceverinol 
Bournaud et al.l I2OIOI: Krumholz fc Thompson! 



2009; 



2012 



20131: IDekel fc Krumholzl 120131) . but it is not clear 
whether they are properly deposited within the disc. 
Evidence against stellar feedback as a direct source 




Figure 16. Inflow within the disc: ratio of gas inflow rate and 
gas mass in the disc, testing eq. I I26II . Representing the left-hand 
side is the median over the simulations (black curve) and the 68 
percentiles (gray shaded area), limited to the snapshots where 
there is net inflow, assumed to undergo VDI. This is compared 
to the right-hand side of the toy model, eq. II26I I. where & and 
arc deduced from the simulations, with v = \/2 and 7 = 1, show- 
ing the median (dashed red curve) and the 68 percentiles (blue 
shaded area). The median of the right-hand side typically under- 
estimates that of the left-hand side by a factor of ~ 2, but the fit 
is acceptable within the scatter and the expected uncertainties, 
which arc especially large at 2 > 3. 
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Figure 17. Inflow within the discs: ratio of gas inflow rate in 
the disc and gas accretion rate into the disc through a sphere 
of 0.1-Rv, testing eq. II27I I. Representing the left-hand side is the 
median over the simulations (black curve) and the 68 percentiles 
(gray shaded area), limited to the snapshots where there is net 
inflow, assumed to undergo VDI. This is compared to the right- 
hand side of the toy model, eq. I I27I I. where 5 is computed from the 
simulations and with 7 = 1, e^f = 0.05, and /^f = 0.5, showing the 
median (dashed red curve) and the 68 percentiles (blue shaded 
area). The fit is good except at 2; > 3, where the uncertainty is 
large. 
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of turbulence is in the observation that the veloc- 
ity dispersion is not tig htly correlated with proximity 
to star-forming clu mps (jForster Schreiber et alj [20091 : 
iGenzel et al.ll20l"l[) . However, outflows can help boost- 
ing up a by simply lowering the surface density of that 
gas tha t is supposed to be stirr ed up by a given energy 
source (jGenel. Dekel fc Cacciat o 2012). Recall that the 
stellar feedback in the simulations used here is under- 
estimated. 

According to eq. ([TS]). the steady-state solution 
implies Mg ~ Mg^ac^sf- Combined with cq. (I26p and 
Tsf = /sfE^^^d from eq. (|20l). this gives in steady state 

A^g,in ^, fsf cr^ ^, fsl ^2 ^27) 

Note that this expression is independent of id- The right 
hand side is of order unity, so the equation implies that 
the rate of gas draining from the disc to the central 
bulge, Mg,in, is expected to be comparable to the rate 
of replenishment by freshly accreting gas, -/Wg,ac- Since 
the steady-state solution, eq. ((T^ . implied that the SFR 
is also comparable to the gas accretion rate, we conclude 
that the inflow rate in the disc is expected to be com- 
parable to the SFR. 

The value of Mg^nAi can serve as a lower limit 
for the growth of mass in the bulge during the period 
At. The ratio of rates in eq. (P7)) can be interpreted as 
an approximation for the bulge-to-total baryonic mass 
ratio in the galaxy. 

The gas mass fraction within the disc radius, Sg, 
can be crudely estimated from the toy models as fol- 
lows. The gas mass is given by the steady-state solu- 
tion, eq. (HH), Mg ~ Mg,ac7"sf ■ Expressing t^i ~ /sfe^f^^d 
and td — 0.007t from eq. (fT4|) . we obtained Tsf/t ~ 0.17. 
We can then write Mg ~ (^sf /t)/g/bMact, where fg is 
the gas fraction in the baryon accretion, /b ~ 0.17 is the 
universal baryon fraction, and il/ac is the total accretion 
rate, such that Mac* — M^. The total mass within the 
disc radius can be crudely estimated as the sum of the 

baryonic mass and the dark mass, Mtot — (/b + A)Mv, 
assuming M(r) oc r within the halo. We thus obtain for 
the gas fraction 

Sg ~ ^^^'/^^^efh < Q ^ (28) 

/b + A 

If only a fraction of the baryons in the halo are within 
the disc radius, the value of dg would be larger accord- 
ingly. When including the cold stars in the disc, the 
value of S could be doubled or higher by a factor of a 
few, so we expect 5 ^ 0.2 — 0.3. 



7.2 Simulations: Disk Inflow 

Figure [TBI and Fig. [T7| test the validity of the toy model 
estimate for the inflow within the disc, eq. (|26l) and 
eq. (|27p respectively. For each equation, we compare 
the left-hand side as measured in the simulations and 
the right-hand side as deduced from the toy model but 



also using quantities that are measured from the simu- 
lations. 

We adopt as our fiducial value for the dissipation 
timescale 7 = 1, and the simulations can actually be 
used to determine the best-fit effective value of 7. 

The cold mass fraction within the disc radius, 5, 
entering the right-hand sides of eq. p6[) and eq. ([27[) . 
is computed for each snapshot inside a sphere of radius 
i?d- "Cold" disk stars are defined by jz/jma^ > 0.7, 
where jz is the specific angular momentum of the star 
particle along the spin axis of the disc, and jmax = 
RVtot is the maximum specific angular momentum at 
the given energy, where R is the radius at the particle 
position and Vtot is the magnitude of the total particle 
velocity. Figure [T8| shows the average and 68% scatter 
of S in the simulations. The median is (5 ~ 0.2 across the 
whole rcdshift range, with a 68% scatter between 0.12 
and 0.25. The value of 6 due to gas only, Sg, is about 
half that value. These measured values are consistent 
with the expectations based on eq. ([28]). 

The dynamical time entering the right-hand-side 
of eq. ([26[) is computed for each snapshot from the sim- 
ulations via eq. ([T4| . Based on eq. ([TT]) . the parameter 

A is derived from the radii as determined from the sim- 
ulations, and V = y/2 is assumed. The Hubble time t is 
approximated by eq. ([T]). 

For the SFR efficiency in eq. ([27|) we use Csf — 0.05. 
This is a high value, estimated to be the effective value 
of egf in our current simulations, dictated by the re- 
quirement that the SFR in the simulations at the given 
resolution roughly matches the Kennicutt relation. This 
effective value of egf is consistent with cq. (IT71) for the 
SFR in the simulations with ts — 4 Myr and where the 
cold gas available for star formation is defined to be the 
mass with temperature lower than lO^K and density 
higher than 1 cm~'^. 

For the analysis here, we have used only snapshots 
where there is a net infiow within the disc, as a signa- 
ture of VDI. Snapshots where the disc radius is smaller 
than 3 kpc are excluded; this minimizes problems re- 
lated to the spatial resolution of the simulations. In- 
cluding these small discs did not make a significant 
difference at z < 3, but did make some difference at 
z = 3 — 4. Since the disc radii at z > 3 are uncertain 
due to wild perturbations, we de-emphasize the results 
in this regime even when the disc radius is larger than 
3 kpc. 

Figure [T6| compares the specific gas infiow rate in 
the disc according to the two sides of eq. ([2^ . The left- 
hand side is measured directly from the simulations, 
while in the right-hand side we make use of the 6 and 
td as computed from the simulations. Overall, the toy 
model is quite successful, at the expected level of un- 
certainty. The average of the left-hand-side is system- 
atically higher than the right-hand side by a factor of 
order two in the range z = 3 — 2. This is at the level 
of accuracy expected from the crude toy model and the 
uncertainties in measuring the simulation results. The 
apparent discrepancy is somewhat larger at z = 4 — 3 
where the uncertainty is larger, and somewhat smaller 
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Figure 18. Steady state in disc: Cold fraction in the disc within 
i?di ^ = A^cold/^tot, Shown is the median over the simulations 
(black curve) and the 68 percentiles (shaded area) for the snap- 
shots where there is net inflow, assumed to undergo VDI. The 
value of (5 ~ 0.2 is consistent with the expectation, eq. II28I I. once 
the mass in cold stars is comparable to the gas mass. 



at z = 2 — 1. A better fit eould be obtained witli either 
a smaller value of 7 or a larger value of v. 

Figure [T71 compares the ratio between the gas inflow 
rate within the disc and the gas accretion rate onto the 
disc as computed by the two sides of eq. ((27|) . The values 
of the parameters used in the right hand side are 7=1, 
esf = 0.05, and /sf = 0.5 The agreement is good except 
in the redshift range z = 3 — 4 where the uncertainty in 
disc radius is large. 



8 STEADY STATE IN THE DISK 

8.1 Toy Model: steady state in disc 

In analogy to the SFR, Afgf = Aig/rgf, based on eq. ([?5|) 
we can express the inflow rate within the disc as Mg in = 
Mg/Tin, where Tin = jia/V)-'^td. The gas-mass conser- 
vation equation, as in eq. (|16p , could then be applied to 
the disc alone, and include draining by three processes, 
namely star formation, outflows, and inflow within the 
disc. Based on eq. (|17l) and eq. (|25|) . each of the cor- 
responding sink terms {i = 1,3) could be expressed in 
the form Mgink.i ~ Mg/ri, where t; = ^j^^td, so the 
continuity equation for the disc is 



g,ac 



AL 



ink : 



sink 



where 



^-1 

'out 



(29) 
(30) 

dH 



and its steady state solution is analogous to eq. 
with Tgf replaced by the smaller timescale r. 

If a steady-state solution exists, we can evaluate the 
steady-state value of (5, the mass fraction of cold disc 
inside i?d, following iDekel. Sari fc Cevering (,2009t ). We 
write 



S = 



(31) 



where /? = Mb/Aftot(^d) is the fraction of baryons 
within the disc radius, estimated to be ~ 0.6. Assum- 
ing that /3 is constant, a time derivative of eq. (|3ip . 
combined with eq. ((29|) . yields 

(32) 

where /g is the mass fraction of gas in the baryonic ac- 
cretion, and where tac = Mb/Mh and isink = Md/Msink- 
Here Afd stands for the cold component in the disc, gas 
and cold young stars. If a steady-state solution exists, 
it will be provided by setting 6 = 0. 

If the sink timescale scales with 6^'^ (see below), 
we write tsink/^ac = bS~^, and obtain from eq. p2p that 
the steady-state solution is the solution of the depressed 
cubic polynomial, 

(33) 



where c = 6/3 /g. Its solution is 



b 

3u ' 



C\l/3 

2 



1/3 



.(34) 



If ^ 1, then c <C 1, and there is no steady-state so- 
lution with (5 ^ 1. If (5 ^ 1, then there is a solution, 
S ~ c/b = /3/g, which is indeed likely to vary rather 
slowly with time. If all three terms of eq. (|33p are com- 
parable, with S ~ 0.3, then we need to appeal to the full 
solution, eq. (|34|) . If the gas dominates the cold disc, we 
can write tsink = Then the long-term time depen- 
dence is via b on a^^. When inserted in eq. p4p . we find 
that the value of 6 is expected to vary rather slowly 
with time, consistent with a steady-state solution. 

The scaling of Tin with is apparent from 

eq. (pS)) . The dependence of Tgt (and therefore Tout) on 
S could be similar if star-formation is limited to the 
disc clumps. Then the SFR is proportional to the mass 
in clumps, namely the product of clump mass and the 
number of clumps, A^cAfc- The Toomre clump m ass is 
Mc/Afd ^ (l/4)(52 (IDekel. Sari fc Ceverinol[2009l) . Con- 
stant values of NcM^, Nc and 6 are consistent with the 
slow variation of 6 implied by eq. (|34|. 



8.2 Simulations: steady state in disc 

Figure [T51 shows that throughout the range z = 3 — 1, 
the average value of d is consistent with being constant, 
as implied by the steady-state solution eq. ([M)) . A slight 
decline may be marginally noticeable toward z ^ 1. The 
average value is 5 ~ 0.2 and the 68% range is roughly 
from 0.12 to 0.24 It is consistent with the crude estimate 
in eq. once the mass in cold stars is comparable to 
the gas mass. 

Figure [TO] shows the evolution of bulge-to-total ra- 
tio. The disc stars are distinguished from bulge stars 
by jz/imax > 0.7. The average is rather constant at 
B/T~ 0.73, with the 68% range from 0.63 to 0.87. This 
is consistent with the steady-state bulge growth by VDI, 
eq. ([?7)) . plus the constancy of S as implied by eq. ([M)) . 



Galaxy- Formation Toy Models 19 



z 



4 3 2 1 




0.6 ~- ^ 
: Simulations at ~ 

: Model, e3,=0.04, 7 = 0.35, (5 = 0.2 E 

Q 5 E \ I \ I I I I \ I I \ I \ \ I I I I I I I I I I I I I I l~ 

0.20 0.25 0.30 0.35 0.40 0.45 0.50 



a 

Figure 19. Steady state in disc: Bulge-to-total ratio within R^. 
Shown are the median (thick black), the average (magenta), and 
the average of the log within the 90% percentiles (green) over 
the simulated galaxies. Shown in comparison is the toy model 
prediction, eq. ^2711, with (5 ~ 0.2, e^f — 0.04 and /^f ~ 0.5. The 
value of B/T and its constancy in time is consistent with the 
steady-state bulge growth by VDI, eq. I I27I I. plus the constancy 
of S as implied by eq. I I34II . 

With S ~ 0.2, Csf ~ 0.04 and fsf ~ 0.5, a value of 
B/T~ 0.7 is predicted for 7 ~ 0.35. This is in the ball- 
park of the value that would yield a good fit in Fig. [T71 
testing the validity of eq. (|26)) for the inflow rate within 
the disc. 



9 CONCLUSION 

We have verified that key processes of galaxy evolution 
in its most active phase, at z = 1 — 4, as captured by 
hydrodynamical simulations in a cosmological context, 
can be approximated by expressions derived from sim- 
ple toy models. The systematic evolution of the stacked 
simulated galaxies is recovered in the best cases at the 
level of 10% accuracy, and in the worse cases only to 
within a factor of two. Even the latter can be useful 
for qualitative understanding and back-of-the-cnvclopc 
estimates. We summarize our main findings below. 

Both the average and the median total specific ac- 
cretion rate into massive haloes are only weakly depen- 
dent on mass and is well approximated as a function of 
redshift by cx (1 -I- z)^/^, eq. ([5]). This includes mergers 
and smoother accretion into the virial radius. The sys- 
tematic growth of virial mass in an individual halo is 
well approximated by a simple exponent in — z, eq. ([9]). 
This implies a slow variation in time of the absolute ac- 
cretion rate into a given halo as it grows in the range 
z = 5—0.3, with a maximum near z 2.2, eq. (jlOp . The 
distribution of specific accretion rate is lognormal, with 
a standard deviation of 0.33 dex. The same expressions 
describe the accretion of baryons in our simulations. 

The penetration of the inflowing baryonic mass {M) 
through the halo into the galaxy inside 0.1 i?v is ^ 50% 



at z = 4 — 2, and somewhat larger at z = 2 — 1, partly 
reflecting the higher stellar fraction in the accretion at 
later times. On the other hand, the specific inflow rate 
M/M at 0.1i?v remains very similar to that at Rv, re- 
flecting the fact that the baryonic mass and its rate of 
change vary in a similar way between the two radii. This 
implies that the toy-model expressions for M/M at i?v 
can serve as a rather accurate estimator of the actual 
characteristic baryon input into the central galaxy. 

The gas in the galaxy obeys a simple mass conserva- 
tion equation, eq. , where the accretion is the source 
and the star formation (and the associated outflows) 
is the si nk. Combined with a universal st ar-formation 
law (e.g. iKrumholz. Dekel fc McKeel[2012[ ). the galaxy 
is predicted to converge to a cosmological quasi-steady 
state, with the SFR following the gas accretion rate into 
the central galaxy and the gas mass constant, eq. p^ . 
We note that this may not be true for less massive 
galaxies and at higher redshifts, where feedback is likely 
to be more effective in suppressin g star formation and 
accumulating a gas rese rvoir (e.g. IZolotov et all 120121 : 
IKrumholz fc DekelllMl) . 

The intense inflows into the dense early galaxies 
trigger violent disc instability (VDI) , which drives mass 
inflow within the disc, building a compact bulge and 
feeding the central black hole. Toy model estimates of 
the inflow rate, performed in several independent ways, 
provide good approximations for the actual inflow rate 
in the simulations, both relative to the gas mass and 
relative to the gas accretion rate, eq. and eq. ([77]) . 
A mass conservation equation that is applied to the disc 
gas alone, eq. ((29|) . predicts a steady state, in which the 
gas in the disc drains by star formation, outflows and 
inflow within the disc to the bulge, and is replenished 
by cosmological accretion. In this steady state, the mass 
fraction of gas and cold disc stars within the disc radius 
is rather constant in time, and the bulge mass is com- 
parable to the disc mass or even higher. In the current 
simulations 5 ~ 0.2 and B/T^ 0.6 — 0.8, while if the gas 
fraction at z 2 is higher, as observed, then 5 could 
be somewhat larger and B/T somewhat smaller (as in 
iDekel. Sari fc Ceverinol 120091) . 

We conclude that these toy models can be very use- 
ful tools in the study of key processes, as seen observa- 
tionally or in simulations, and can thus help us gain 
understanding of the complex process of galaxy forma- 
tion. In turn, the toy models allow us to perform simple 
analytic calculations to predict the expected evolution 
of galaxies. A more detailed treatment of star forma- 
tion, and the effects of very strong outflows, are being 
and will be incorporated into the simulations, and are 
yet to be reflected in further testing of toy models. 
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APPENDIX A: USEFUL RELATIONS 

We summarize here the more accurate cosmological 
relations, valid in the standard ACDM cosmology 
at all redshifts, which the toy-model expressions 
of this paper approximate in the EdS regime, 
z > 1. This is largely basic material, b as ed fo r 
example on Lah av. Rees. Li lje. fc Priraack (11991 ): 
iCarroU, Press. fcTurner (1993) and iMofcWhiti 



('2002 
(j200e 
from 



Birnboim fc Dekell (|'200l : iDekel fc BirnboimI 
with the addition of accre tion rate s 
Neistein. van den Bosch fc Dekell ()2006D : 



INeistein fc Dekell (|2008bf r 



Al Cosmology 

The basic parameters characterizing a flat cosmological 
model in the matter era are the current values of the 
mean mass density parameter fim and the Hubble con- 
stant Hq . At the time associated with expansion factor 
a = l/(l + z), the vacuum-energy density parameter is 
17a (a) = 1 ~ 17m(a) and 

The Hubble constant is 

Hia) = Ho{nA + Q^a-y^^, (A2) 

and the age of the universe is 

. 2 , sinh'^(|l-a„(a)|/a„(a))i/^ 

'^"^^S"^^"^ (|l-an(a)|)V2 ■ ^"^'^ 

The mean mass density is 

Pu ~ 2.76 X 10-3°a„o.3 ^0.7 , (A4) 

where il,„Q 3 = f2ni/0.3, h = _ffo/100kms^^ Mpc^^, and 
ho.7 = h/o'.T. 

A2 Virial relations 

The virial relations between halo mass, velocity and ra- 
dius, 



GNU 



(47r/3)i?3 



(A5) 



i?ioo ^ 1.03Afif ^1/3 , (A6) 



become 

^200 - imMli^A-jl'' 

where M12 = M^/W^Mq, = K/200kms-i, 

i?ioo = -Rv/lOOkpc, and 

A={^2m^^o.zhlj)-^''^a. (A7) 

An approximation for A (a) in a flat universe (Bryan & 
Norman 1998) is: 

A(a) ~ (187r2 - 82r2A(a) - 3917A(a)^)/an(a) ■ (A8) 

In the EdS regime, z > 1, or when referring to i?2oo 
instead of at all redshifts, A ~ (1 -\- z)~^ to an accu- 
racy of a few percent. The virial temperature, defined by 
kTy/m = (1/2)V^^. For an isotropic, isothermal sphere. 
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this equals cr^, where a is the one-dimensional veloc- 
ity dispersion and the internal energy per unit mass is 
e = (3/2)cr2. Thus 



Tg ~ 2.87X^200 ' 
where Tq = Tv/IO^K. 

A3 Press Schechter 



(A9) 



Linear fluctuation growth is given by (iLahav et al.l 
119911: ICarroll. Press fc Turnedll992t[Mo fc Whitell200l 



D{a) = 
where 
9{a) ' 



.9(a) 
5(1 



(AlO) 



(All) 



n, 



o (n) I (l + "m(a)/2) 

(1 + i2A(a)/70) 



The CDM pow er spectrum is crudely approximated 
by CBardeen et al.i,1986.f l 

(A12) 
(A13) 

(5.46(7)3 + (6.71g)*]-l/^ 
(A14) 



P{k) cx fcT2(fc), 
with 

^ ln(l + 2.34,) 
^ ^ 2.34 9 

X [1 -h3.89g+ (16.1g)2 

where 

q = k/{n„,h'^Mpc-^). 
It is normalized by erg at i? = 8/i^^Mpc, where 

1 r°° 

o-2(i?) = — / dkk^ P{k)W^{kR), 
27r ./n 



(A15) 



and with the Fourier transform of the top-hat window 
function 



W{x) — 3(sina; — xcosx)/x^ 



(A16) 



In the Press Schechter (PS) approximation, the 
characteristic halo mass Afps(a) is defined as the mass 
of the l-a fluctuation, 



1 = iy{M, a) 



1.68, 



(A17) 



D{a) a{M) 

where M and the comoving radius R are related via the 
universal density today: M = ^Po R^- The mass of 2- 
cr fluctuations is obtained by setting t/(M, a) = 2, etc. 
Based on the improved formalism of ISheth fc TormenI 
()2002|) . the fraction of total mass in haloes of masses 
exceeding M is 

^^(>M,a)^ 0.4 (1 + MVrfe ^0-85- 
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Figure Al. The Prcss-Schcchtcr mass as a function of redshift 
(blue). Also shown are the approximations log A/ps = 12.7 — 1.29z 
(green), and logMps = 12.78 - 1.46^°-** (red). 



This fraction for l-a, 2-a, and 3-(T fluctuations is 22%, 
4.7%, and 0.54% respectivelyEI 

Figure lAll shows the PS mass Mpg as a function of 
redshift. For the standard ACDM with f2,„ = 0.27 and 
as = 0.8 its value at z = is A/ps,o ^ 6 x lO^^M©. 
A practical fit for z ^ 4 with ~ 10% accuracy is pro- 
vided by a power law in this semi-log plot: logMpg = 
12.7— 1.292. At larger redshifts this gradually becomes 
an underestimate. A fit that is accurate to 20% in z < 6 
and is an underestimate by a factor of 2 at z = 8 is 
log A/ps = 12.78 - 1.46z°-^^. Trying to provide crude 
power-law approximations, we find that A/ps oc a^'^ oc 
^3.5 g^j.g Qj-mJe approximations in the range 5j z ^ 1, 
and that A/ps oc a"^ oc are good to within a factor 
of 2 in the range ^ z ^ 2. These power laws become 
overestimates at higher redshifts. 



A4 Accretion Rate 

The average growth of a halo, following its main 
progenitor, has been derived by a fit to merger 
trees from the Millennium c osmological simula- 
tion (jSpringel fc HernquistI 120051 ). using a functional 
form that is motivated by the EPS approximation 
(jNeistein fc Dekel l2008bl ). In terms of the self-similar 
dimensionless time variable lu = Sc[D{a)~^ — 1], the 
growth rate of halo of mass A/ at ui is 

dM 
duj 

and the mass at u, given that the mass at z = is A/g, 
is 



(A19) 



/ 0.4 \ ^ /0.85i^\ 1^ 

[' + ^ j [-Vrj ■ ^^^'^ MHMo) = (A/-^ + a/3c.)- V. , 



(A20) 



^ Note that in the initial conditions for our simulations we use 
the more accurate CAMB numerical solutions based on CAMfast 
(Scljak and Zaldarriaga, |http: / /lambda.gsfc.nasa.gov/toolbox/[ l. 



^ Note that the Sheth-Tormen approximation be- 
comes quite inaccurate at ve ry high redshifts 
l|KlvDin. Truiillo-Gomez fc Primac^l2011^ . 
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where M is in lO^^Mg. The value /3 = 0.14 fits well the 
Millennium trees for M in the range 10^^ — IO^'^A/q, 
despite the fact that according to EPS it should vary 
as ^ = (n + 3)/6, where n is the local power-spectrum 
index. At lower masses, /? should be smaller, so the mass 
dependence of the accretion rate becomes even closer 
to linear. For erg = 0.81, the best fit is with a = 0.623 
(it scales as oc u^^). The accuracy of this fit to the 
Millennium trees is better than 5% for z < 5. 

The following approximations (better than 0.5% ac- 
curacy at all z for 51„i = 0.27, h = 0.7, o-g — 0.81) relate 
uj to z, 



1.28[(1 



0.086(1 + 



0.22e" 



(A21) 



and provide the standard accretion rate M = udM / duj 
with 

UJ = -0.0476[(l + z)+0.093(l-fz)-i-22]2-5Gyr-i.(A22) 

The latter scales in proportion to oc /io.7- Note that in 
the EdS regime, valid at high z, the asymptotic depen- 
dence is w cx (l-l-z)^-^, because D{t) = (l-|-z)~^ oc i^/^, 
so M cx t-^/'^ cx (1 + z)5/2. 

Based on these expressions, a practical approxima- 
tion for the typical baryonic accretion rate for a halo of 
mass M at z is 



Mb = 80 (1 + z)^ /o.i7 Mq yr- 



(A23) 



where /0.17 is the baryonic fraction in units of 0.17. 
With the EdS asymptotic value ji = 2.5, the accuracy 
is < 5% for z > 1, and it becomes an underestimate of 
~ 20% at z = 0. With /i = 2.4 the accuracy is < 5% 
for 0.2 < z < 5, and is an underestimate of ^ 10% at 
z = and z = 10. For most practical purposes we adopt 
/3 = 0.14 and = 2.4. 

Figure IA2I is a useful summary for appreciating the 
average baryon accretion rate onto haloes as they grow 
during cosmological evolution. Each of the blue curves 
refers to a given value of Mb (as marked). They are 
obtained using the inverse of eq. ([Hoi) and cq. ((X22)) . 
Each of the red curves shows the average mass growth of 
the main progenitor of a halo that ends up with a given 
mass Mo at z = (for logMo = 11.5, 12.5, 13.5), using 
eq. (|A20p and eq. (jA2ip . For a given halo of Mv at z, one 
can read the baryon accretion rate at the correspond- 
ing point in the Mv , z plane by its position relative to 
the nearest reference blue curve. For example, a halo 
of Mv = 10^2 Mq at z = 2 accretes baryons at about 
8OM0yr-\ as in eq. 0. 

By following the red curves one can read the history 
of accretion rate onto a given halo as it grows. Take for 
example the middle red curve, referring to a halo of 
Mo = 3 X IO^^Mq at z = 0. It grew between z = 4 
and 0.5 from Mv ~ 2 x IO^Mq to 2 x IO^^^q, while 
maintaining an accretion rate of 40 — 60 Mq yr~^. 

A5 Clump Migration by Dynamical Friction 

The timescale for inward clump migration in a 
m arginally unstable disc, Q ~ 1, has been estimated 
bv lDekel. Sari fc CeveriMi (|2009f ) to roughly be 



> 
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Figure A2. Baryonic accretion rate, halo mass and rcdshift. The 
sohd blue curves refer to constant values of baryon accretion rate 
(an upper limit to the SFR) as marked. The long-dashed red 
curves show the average mass growth of the main progenitor of 
a halo of a given mass Mq at 2; = 0, for logMg = 11.5, 12.5, 13.5 
from bottom to top. The short-dash green curve is the Press- 
Schechter mass. 



(A24) 

where S is the cold mass fraction defined by Md = 
(SMtot, with Md the cold mass in the disc and Mtot 
the total mass encompassed by the disc radius i?d- The 
timescale for evacuating a disc mass Md is then 

^cvac — ^mig 5 (A25) 

where rric ^ 0.2 is the instantaneous fraction of the disc 
mass in clumps. This estimate was based on the en- 
ergy exchange in the gravitational encounters between 
clumps. 

Here we provide an alternative crude estimate of 
the clump migration time due to the dynamical fric- 
tion that is exerted on the clump by the off-clump disc 
mass. Consider a clump of mass Mc in a circular or- 
bit with velocity V ~ {GMtot/ Rd)^^'^ ■ The deceleration 
due to dynamical fr iction is approximated by t he Chan- 
drasekhar formula (jBinnev fc Tremainell2008l sec. 8.1) 
to be 



V c^-AirG^ Ah PA (InA)y- 



(A26) 



In a Q 1 marginally unstable disc, the clump mass is 
related to the disc mass by 

Me ~ (J^Md = 53 Mtot , (A27) 

and the disc half thickness h is related to its radius by 

h/RA - 5 (A28) 

(e.g. iDekel. Sari fc Ceverinol I2OO90 . Thus, the mean 
density in the disc is 

Md Mtot 



Pd 



2TTRlh 



In the Coulomb logarithm A ~ M^/Mc 
assumed to be ~ 2. We finally obtain 



(A29) 
5^^ so In A is 
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imig = V/V ^- 6-hd . (A30) 

This is indeed comparable to the estimate in eq. (|A24p 
for 6 ~ 0.2 — 0.3, the typical val ue in the s teady-state 
solution of high-z disks (IPekeh Sari &: Ceve rino 2Q0l; 
ICacciato. Dekel fc Genelll2012h and in our current sim- 
ulations, Fig. [m The disc evacuation time by clump 
migration due to dynamical friction is thus 

^cvac - S-hd . (A31) 



